
#include <thrust/async/for_each.h>
#include <thrust/iterator/counting_iterator.h>

#include <algorithm/analytical_model/Constants.h>
#include <algorithm/analytical_model/algorithm/cpu/ConvolutionMethod.h>
#include <algorithm/analytical_model/algorithm/cuda/ConvolutionMethod.cuh>
#include <algorithm/analytical_model/algorithm/cuda/GeometryUtils.cuh>
#include <algorithm/analytical_model/algorithm/cuda/Utils.cuh>
#include <algorithm/analytical_model/algorithm/cuda/Utils.h>
#include <algorithm/analytical_model/algorithm/cuda/VectorMath.cuh>
#include <algorithm/analytical_model/algorithm/cuda/inter_op/Data.h>
#include <algorithm/analytical_model/algorithm/cuda/math/fast_math/Erff.cuh>
#include <algorithm/analytical_model/algorithm/dnn/params/UNIZAR2.h>
#include <algorithm/analytical_model/algorithm/precomp/GlParamArray.h>
#include <utils/Format.h>
#include <utils/Utils.h>
#include <utils/config.h>
#include <utils/cuda/containers/Vector.cuh>

#include <cuda_runtime.h>
#include <cuda_runtime_api.h>

#include <algorithm>
#include <bit>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <iostream>
#include <memory>
#include <numbers>
#include <stdexcept>

using namespace solar::cuda;
using namespace solar::cuda::inter_op;

namespace
{
    __constant__ constexpr float cm_gl_param_10_precomp[] = {
        0.00000F,  2.00000F,  -0.57735F, 0.57735F,  1.00000F,  1.00000F,  -0.77460F, 0.00000F,
        0.77460F,  0.55556F,  0.88889F,  0.55556F,  -0.86114F, -0.33998F, 0.33998F,  0.86114F,
        0.34785F,  0.65215F,  0.65215F,  0.34785F,  -0.90618F, -0.53847F, 0.00000F,  0.53847F,
        0.90618F,  0.23693F,  0.47863F,  0.56889F,  0.47863F,  0.23693F,  -0.93247F, -0.66121F,
        -0.23862F, 0.23862F,  0.66121F,  0.93247F,  0.17132F,  0.36076F,  0.46791F,  0.46791F,
        0.36076F,  0.17132F,  -0.94911F, -0.74153F, -0.40585F, 0.00000F,  0.40585F,  0.74153F,
        0.94911F,  0.12948F,  0.27971F,  0.38183F,  0.41796F,  0.38183F,  0.27971F,  0.12948F,
        -0.96029F, -0.79667F, -0.52553F, -0.18343F, 0.18343F,  0.52553F,  0.79667F,  0.96029F,
        0.10123F,  0.22238F,  0.31371F,  0.36268F,  0.36268F,  0.31371F,  0.22238F,  0.10123F,
        -0.96816F, -0.83603F, -0.61337F, -0.32425F, 0.00000F,  0.32425F,  0.61337F,  0.83603F,
        0.96816F,  0.08127F,  0.18065F,  0.26061F,  0.31235F,  0.33024F,  0.31235F,  0.26061F,
        0.18065F,  0.08127F,  -0.97391F, -0.86506F, -0.67941F, -0.43340F, -0.14887F, 0.14887F,
        0.43340F,  0.67941F,  0.86506F,  0.97391F,  0.06667F,  0.14945F,  0.21909F,  0.26927F,
        0.29552F,  0.29552F,  0.26927F,  0.21909F,  0.14945F,  0.06667F};
    ;

    /**
     * 高斯勒让德数值积分实现，默认针对矩形形进行计算
     * 将[-delta_x,delta_x]x[-delta_y,delta_y]的区域映射到[-1,1]x[-1,1]上
     *
     * @param m 在x方向上取点个数
     * @param n 在y方向上取点个数
     * @param x,y 被积分点位置坐标
     * @param xx,yy 待积区域左上角坐标
     * @param delta_x, delta_y 小平行四边形边长的一半
     * @param sigma1, sigma2 UNIZAR/HUANG/CAUCHY中模型的参数
     * @param mod 对应三种模型，前两者公式相同
     * @param flux_param 卷积模型中单点反射能量总值
     * @param gl_params 10行，每一行i*2个浮点数，前i个为节点坐标，后i个为对应系数
     * @return
     */
    inline __device__ auto gaussLegendre(int m, int n, float x, float y, float xx, float yy,
                                         float delta_x, float delta_y, float sigma1, float sigma2,
                                         solar::ConvolutionModelType mod, float flux_param
                                         /*float* gl_params*/) -> float
    {
        const float* gl_params = cm_gl_param_10_precomp;
        int gl_p1 = m * (m - 1);
        int gl_p2 = n * (n - 1);
        float sum = 0;
        for (int i = 0; i < m; i++)
        {
            float tmp_sum = 0;
            float xi = gl_params[gl_p1 + i];
            float ai = gl_params[gl_p1 + m + i];
            for (int j = 0; j < n; j++)
            {
                float xj = gl_params[gl_p2 + j];
                float aj = gl_params[gl_p2 + n + j];
                float u = xx + xi * delta_x - x;
                float v = yy + xj * delta_y - y;
                // 暂时使用UNIZAR的参数
                //            tmp_sum += aj * delta_y;
                if (mod == solar::ConvolutionModelType::kUNIZAR)
                {
                    tmp_sum += aj * delta_y * 0.5F / inter_op::kMathPi / sigma1 / sigma2 *
                               __expf(-0.5F * (u * u / sigma1 / sigma1 + v * v / sigma2 / sigma2));
                }
                else
                {
                    tmp_sum += aj * delta_y * sigma1 * sigma1 / inter_op::kMathPi /
                               (1 + sigma1 * sigma1 * (u * u + v * v)) /
                               (1 + sigma1 * sigma1 * (u * u + v * v));
                }
            }
            sum += tmp_sum * ai * delta_x;
        }
        return sum * flux_param;
        //    return sum;
    }

    /**
     * 高斯勒让德数值积分实现，默认针对平行四边形形进行计算
     * 将[-delta_x,delta_x]x[-delta_y,delta_y]的区域映射到[-1,1]x[-1,1]上
     *
     * @param m 在x方向上取点个数
     * @param n 在y方向上取点个数
     * @param x,y 被积分点位置坐标
     * @param xx,yy 待积区域左上角坐标
     * @param delta_x, delta_y 小平行四边形边长的一半
     * @param sigma1, sigma2 UNIZAR/HUANG/CAUCHY中模型的参数
     * @param mod 对应三种模型，前两者公式相同
     * @param flux_param 卷积模型中单点反射能量总值
     * @param gl_params 10行，每一行i*2个浮点数，前i个为节点坐标，后i个为对应系数
     * @return
     */
    inline __device__ auto gaussLegendre(int m, int n, float x, float y, float xx, float yy,
                                         float2 delta_x, float2 delta_y, float sigma1, float sigma2,
                                         solar::ConvolutionModelType mod, float flux_param
                                         /*float* gl_params*/) -> float
    {
        const float* gl_params = cm_gl_param_10_precomp;
        int gl_p1 = m * (m - 1);
        int gl_p2 = n * (n - 1);
        float sum = 0;
        for (int i = 0; i < m; i++)
        {
            float tmp_sum = 0;
            float xi = gl_params[gl_p1 + i];
            float ai = gl_params[gl_p1 + m + i];
            for (int j = 0; j < n; j++)
            {
                float xj = gl_params[gl_p2 + j];
                float aj = gl_params[gl_p2 + n + j];
                /* 需要考虑在两个方向上的移动对x和y的影响 */
                float u = xx + (xi * delta_x).x + (xj * delta_y).x - x;
                float v = yy + (xj * delta_y).y + (xi * delta_x).y - y;
                // 暂时使用UNIZAR的参数
                //            tmp_sum += aj * delta_y;
                // Cauchy
                //            tmp_sum += aj * sigma1 * sigma1 / MATH_PI /
                //            pow(1+sigma1*sigma1*(u*u+v*v),2);
                // UNIZAR
                if (mod == solar::ConvolutionModelType::kiCauchy)
                {
                    tmp_sum += aj * sigma1 * sigma1 / inter_op::kMathPi /
                               __powf(1 + sigma1 * sigma1 * (u * u + v * v), 2);
                }
                else
                {
                    tmp_sum += aj * 0.5F / inter_op::kMathPi / sigma1 / sigma2 *
                               __expf(-0.5F * (u * u / sigma1 / sigma1 + v * v / sigma2 / sigma2));
                }
                //            if (mod == UNIZAR)
                //                tmp_sum += aj * 0.5f / MATH_PI / sigma1 / sigma2
                //                           * exp(-0.5f* (u * u / sigma1 / sigma1 + v * v/ sigma2/
                //                           sigma2));
                //            else {
                //                tmp_sum += aj * sigma1 * sigma1 / MATH_PI /
                //                pow(1+sigma1*sigma1*(u*u+v*v),2);
                //            }
            }
            // abs函数中的是Jacobian行列式
            sum += tmp_sum * ai * std::abs(delta_y.x * delta_x.y - delta_y.y * delta_x.x);
        }
        return sum * flux_param;
        //    return sum;
    }

    __global__ void convCylinder(float* d_array, float theta_0, float sigma1, float sigma2,
                                 float flux_param, int2 recv_row_col, float3 recv_size,
                                 float3 recv_pos, int2 helio_row_col, float3 helio_vert0,
                                 float3 delta_x, float3 delta_y, float3 focus_center, float3 x_axis,
                                 float3 y_axis, float3 z_axis, float* gl_params, bool* d_micro_sb,
                                 solar::ConvolutionModelType mode)
    {
        std::int64_t my_id = getThreadID();
        //    int2 row_col = hflcal->getRowCol();
        if (my_id >= recv_row_col.x * recv_row_col.y)
        {
            return;
        }
        // 计算当前点对应接收器平面中网格点坐标i,j
        int2 pos_ij = make_int2(my_id % recv_row_col.x, my_id / recv_row_col.x);
        // 计算网格点的世界坐标
        float pixel_angle =
            theta_0 + 1.0F * (pos_ij.x + 0.5F) / recv_row_col.x * 2 * inter_op::kMathPi;

        const auto project_cos =
            dot(z_axis, make_float3(__cosf(pixel_angle), 0, __sinf(pixel_angle)));
        if (project_cos <= 0)
        {
            return;
        }

        float3 pos =
            recv_pos +
            make_float3(__cosf(pixel_angle) * recv_size.x, 0, __sinf(pixel_angle) * recv_size.x) +
            make_float3(0, 1.0F * (pos_ij.y + 0.5F) / recv_row_col.y * recv_size.y, 0);
        // 计算网格点在成像平面投影的世界坐标
        float3 projected_pos;
        // 斜平行投影找对应点
        std::tie(projected_pos, std::ignore) = intersect(z_axis, focus_center, pos, -z_axis);
        // 投影点在成像平面上的二维坐标
        float u = dot(x_axis, projected_pos - focus_center);
        float v = dot(y_axis, projected_pos - focus_center);

        float sum = 0;

        for (int i = 0; i < helio_row_col.y; i++)
        {
            for (int j = 0; j < helio_row_col.x; j++)
            {
                if (d_micro_sb[i * helio_row_col.x + j])
                {
                    continue;
                }
                float3 left_up = helio_vert0 + delta_x * j / helio_row_col.x +
                                 delta_y * i / helio_row_col.y - focus_center;
                float3 right_bottom = helio_vert0 + delta_x * (j + 1) / helio_row_col.x +
                                      delta_y * (i + 1) / helio_row_col.y - focus_center;
                float3 left_bottom = helio_vert0 + delta_x * j / helio_row_col.x +
                                     delta_y * (i + 1) / helio_row_col.y - focus_center;
                float xx = dot(left_up + right_bottom, x_axis) * 0.5F;
                float yy = dot(left_up + right_bottom, y_axis) * 0.5F;
                float3 tmp_x = left_bottom - left_up;
                float3 tmp_y = right_bottom - left_bottom;
                float2 delta_xx = make_float2(dot(tmp_x, x_axis), dot(tmp_x, y_axis)) * 0.5F;
                float2 delta_yy = make_float2(dot(tmp_y, x_axis), dot(tmp_y, y_axis)) * 0.5F;
                //            float delta_xx = dot(right_bottom-left_up, x_axis) * 0.5;
                //            float delta_yy = dot(right_bottom-left_up, y_axis) * 0.5;
                //            sum +=
                //            gauss_legendre(6,6,u,v,xx,yy,delta_xx,delta_yy,sigma1,sigma2,UNIZAR,flux_param,gl_params);
                //            sum +=
                //            gauss_legendre(6,6,u,v,xx,yy,delta_xx,delta_yy,sigma1,sigma2,Cauchy,flux_param,gl_params);
                sum += gaussLegendre(4, 4, u, v, xx, yy, delta_xx, delta_yy, sigma1, sigma2, mode,
                                     flux_param /*, gl_params*/);
            }
            // sum += tmp_sum;
        }

        sum *= project_cos;
        atomicAdd(&(d_array[my_id]), sum);
    }

    __global__ void convCylinderFast(float* d_array, float theta_0, float sigma1, float sigma2,
                                     float flux_param, int2 recv_row_col, float3 recv_size,
                                     float3 recv_pos, int2 helio_row_col, float3 focus_center,
                                     float3 x_axis, float3 y_axis, float3 z_axis, bool* d_micro_sb,
                                     float2* d_new_rectangle, float jacobian_coeff,
                                     float gradient_inv, solar::ConvolutionModelType mode)
    {
        const std::int64_t my_id = getThreadID();
        //    int2 row_col = hflcal->getRowCol();
        if (my_id >= recv_row_col.x * recv_row_col.y)
        {
            return;
        }
        // printf("%f %f %f\n", x_axis.x, x_axis.y, x_axis.z);
        // printf("%f   %f\n", sigma1, sigma2);
        // 计算当前点对应接收器平面中网格点坐标i,j
        const int2 pos_ij = make_int2(my_id % recv_row_col.x, my_id / recv_row_col.x);
        // 计算网格点的世界坐标
        const float pixel_angle =
            theta_0 + 1.0F * (pos_ij.x + 0.5F) / recv_row_col.x * 2.0F * inter_op::kMathPi;

        const auto project_cos =
            dot(z_axis, make_float3(__cosf(pixel_angle), 0, __sinf(pixel_angle)));
        if (project_cos <= 0)
        {
            return;
        }

        const float3 pos =
            recv_pos +
            make_float3(__cosf(pixel_angle) * recv_size.x, 0, __sinf(pixel_angle) * recv_size.x) +
            make_float3(0, 1.0F * (pos_ij.y + 0.5F) / recv_row_col.y * recv_size.y, 0);
        // 计算网格点在成像平面投影的世界坐标
        float3 projected_pos;
        // 斜平行投影找对应点
        std::tie(projected_pos, std::ignore) = intersect(z_axis, focus_center, pos, -z_axis);
        // 投影点在成像平面上的二维坐标
        float u = dot(x_axis, projected_pos - focus_center); // x0
        float v = dot(y_axis, projected_pos - focus_center); // y0

        const float x13 = d_new_rectangle[2].x;
        const float x10 = d_new_rectangle[1].x;
        const float y11 = d_new_rectangle[1].y;
        const float y10 = d_new_rectangle[0].y;
        const float sigma1_denominator = 1.0F / (std::numbers::sqrt2_v<float> * sigma1);
        const float sigma2_denominator = 1.0F / (std::numbers::sqrt2_v<float> * sigma2);
        const auto actual_u = u + gradient_inv * (v - y11);
        float sum = std::abs(1.0F * 0.25F *
                             (fast_math::erff((x13 - actual_u) * sigma1_denominator) -
                              fast_math::erff((x10 - actual_u) * sigma1_denominator)) *
                             (fast_math::erff((y11 - v) * sigma2_denominator) -
                              fast_math::erff((y10 - v) * sigma2_denominator)) *
                             jacobian_coeff * flux_param);
        // printf("sum=%f jacobian_coeff=%f flux_param=%f\n", sum, jacobian_coeff, flux_param);
        // printf("%f, %f, %f, %f, %f, %f, %f, %f, %f, %f\n", sigma1, sigma2, flux_param, x13, x10,
        // y11, y10, u, v, sum);

        const auto delta_x = std::abs(x13 - x10);
        const auto delta_y = std::abs(y10 - y11);

        for (int i = 0; i < helio_row_col.y; i++)
        {
            for (int j = 0; j < helio_row_col.x; j++)
            {
                if (!d_micro_sb[i * helio_row_col.x + j])
                {
                    continue;
                }
                const auto x13_1 = x10 + delta_x * ((j + 1) / helio_row_col.x);
                const auto x10_1 = x10 + delta_x * (j / helio_row_col.x);
                const auto y11_1 = y10 + delta_y * ((i + 1) / helio_row_col.y);
                const auto y10_1 = y10 + delta_y * (i / helio_row_col.y);
                const auto actual_u_1 = u + gradient_inv * (v - y11);
                sum -= std::abs(1.0F * 0.25F *
                                (fast_math::erff((x13_1 - actual_u_1) * sigma1_denominator) -
                                 fast_math::erff((x10_1 - actual_u_1) * sigma1_denominator)) *
                                (fast_math::erff((y11_1 - v) * sigma2_denominator) -
                                 fast_math::erff((y10_1 - v) * sigma2_denominator)) *
                                1.0F * flux_param);
            }
        }
        // printf("%d\n", sum);

        sum *= project_cos;
        atomicAdd(&(d_array[my_id]), sum);
    }

    __global__ void convCylinderFast2(float* d_array, float theta_0,
                                      dnn::UNIZAR2Output model_output, int2 recv_row_col,
                                      float3 recv_size, float3 recv_pos, int2 helio_row_col,
                                      float3 focus_center, float3 x_axis, float3 y_axis,
                                      float3 z_axis, bool* d_micro_sb, float jacobian_coeff,
                                      float gradient_inv, solar::ConvolutionModelType mode)
    {
        const std::int64_t my_id = getThreadID();
        //    int2 row_col = hflcal->getRowCol();
        if (my_id >= recv_row_col.x * recv_row_col.y)
        {
            return;
        }
        // printf("%f %f %f\n", x_axis.x, x_axis.y, x_axis.z);
        // printf("%f   %f\n", sigma1, sigma2);
        // 计算当前点对应接收器平面中网格点坐标i,j
        const int2 pos_ij = make_int2(my_id % recv_row_col.x, my_id / recv_row_col.x);
        // 计算网格点的世界坐标
        const float pixel_angle =
            theta_0 + 1.0F * (pos_ij.x + 0.5F) / recv_row_col.x * 2.0F * inter_op::kMathPi;

        const auto project_cos =
            dot(z_axis, make_float3(__cosf(pixel_angle), 0, __sinf(pixel_angle)));
        if (project_cos <= 0)
        {
            return;
        }

        const float3 pos =
            recv_pos +
            make_float3(__cosf(pixel_angle) * recv_size.x, 0, __sinf(pixel_angle) * recv_size.x) +
            make_float3(0, 1.0F * (pos_ij.y + 0.5F) / recv_row_col.y * recv_size.y, 0);
        // 计算网格点在成像平面投影的世界坐标
        float3 projected_pos;
        // 斜平行投影找对应点
        std::tie(projected_pos, std::ignore) = intersect(z_axis, focus_center, pos, -z_axis);
        // 投影点在成像平面上的二维坐标
        const float u = dot(x_axis, projected_pos - focus_center); // x0
        const float v = dot(y_axis, projected_pos - focus_center); // y0

        const float x13 = model_output.x3_;
        const float x10 = model_output.x0_;
        const float y11 = model_output.y1_;
        const float y10 = model_output.y0_;
        const float sigma1_denominator =
            1.0F / (std::numbers::sqrt2_v<float> * model_output.sigma_x_);
        const float sigma2_denominator =
            1.0F / (std::numbers::sqrt2_v<float> * model_output.sigma_y_);
        const auto actual_u = u + gradient_inv * (v - y11);
        float sum = std::abs(1.0F * 0.25F *
                             (fast_math::erff((x13 - actual_u) * sigma1_denominator) -
                              fast_math::erff((x10 - actual_u) * sigma1_denominator)) *
                             (fast_math::erff((y11 - v) * sigma2_denominator) -
                              fast_math::erff((y10 - v) * sigma2_denominator)) *
                             1.0F * model_output.peak_flux_param_);

        const auto delta_x = std::abs(x13 - x10);
        const auto delta_y = std::abs(y10 - y11);

        for (int i = 0; i < helio_row_col.y; i++)
        {
            for (int j = 0; j < helio_row_col.x; j++)
            {
                if (!d_micro_sb[i * helio_row_col.x + j])
                {
                    continue;
                }
                const auto x13_1 = x10 + delta_x * ((j + 1) / helio_row_col.x);
                const auto x10_1 = x10 + delta_x * (j / helio_row_col.x);
                const auto y11_1 = y10 + delta_y * ((i + 1) / helio_row_col.y);
                const auto y10_1 = y10 + delta_y * (i / helio_row_col.y);
                const auto actual_u_1 = u + gradient_inv * (v - y11);
                sum -= std::abs(1.0F * 0.25F *
                                (fast_math::erff((x13_1 - actual_u_1) * sigma1_denominator) -
                                 fast_math::erff((x10_1 - actual_u_1) * sigma1_denominator)) *
                                (fast_math::erff((y11_1 - v) * sigma2_denominator) -
                                 fast_math::erff((y10_1 - v) * sigma2_denominator)) *
                                1.0F * model_output.peak_flux_param_);
            }
        }
        // printf("%d\n", sum);

        sum *= project_cos;
        atomicAdd(&(d_array[my_id]), sum);
    }

    __global__ void convCylinderFastWithMerge(float* d_array, float theta_0, float sigma1,
                                              float sigma2, float flux_param, int2 recv_row_col,
                                              float3 recv_size, float3 recv_pos, int2 helio_row_col,
                                              float3 helio_vert0, float3 delta_x, float3 delta_y,
                                              float3 focus_center, float3 x_axis, float3 y_axis,
                                              float3 z_axis, float* gl_params, bool* d_micro_sb,
                                              float2* d_new_rectangle, float jacobian_coeff,
                                              solar::ConvolutionModelType mode)
    {
        const std::int64_t my_id = getThreadID();
        //    int2 row_col = hflcal->getRowCol();
        if (my_id >= recv_row_col.x * recv_row_col.y)
        {
            return;
        }
        // printf("%f   %f\n", sigma1, sigma2);
        // 计算当前点对应接收器平面中网格点坐标i,j
        const int2 pos_ij = make_int2(my_id % recv_row_col.x, my_id / recv_row_col.x);
        // 计算网格点的世界坐标
        const float pixel_angle =
            theta_0 + 1.0F * (pos_ij.x + 0.5F) / recv_row_col.x * 2.0F * inter_op::kMathPi;

        const auto project_cos =
            dot(z_axis, make_float3(__cosf(pixel_angle), 0, __sinf(pixel_angle)));
        if (project_cos <= 0)
        {
            return;
        }

        const float3 pos =
            recv_pos +
            make_float3(__cosf(pixel_angle) * recv_size.x, 0, __sinf(pixel_angle) * recv_size.x) +
            make_float3(0, 1.0F * (pos_ij.y + 0.5F) / recv_row_col.y * recv_size.y, 0);
        // 计算网格点在成像平面投影的世界坐标
        float3 projected_pos;
        // 斜平行投影找对应点
        std::tie(projected_pos, std::ignore) = intersect(z_axis, focus_center, pos, -z_axis);
        // 投影点在成像平面上的二维坐标
        float u = dot(x_axis, projected_pos - focus_center); // x0
        float v = dot(y_axis, projected_pos - focus_center); // y0

        float x13 = d_new_rectangle[2].x;
        float x10 = d_new_rectangle[1].x;
        float y11 = d_new_rectangle[1].y;
        float y10 = d_new_rectangle[0].y;
        const float sigma1_denominator = 1.0F / (std::numbers::sqrt2_v<float> * sigma1);
        const float sigma2_denominator = 1.0F / (std::numbers::sqrt2_v<float> * sigma2);
        float sum = std::abs(1.0F * 0.25F *
                             (fast_math::erff((x13 - u) * sigma1_denominator) -
                              fast_math::erff((x10 - u) * sigma1_denominator)) *
                             (fast_math::erff((y11 - v) * sigma2_denominator) -
                              fast_math::erff((y10 - v) * sigma2_denominator)) *
                             jacobian_coeff * flux_param);
        // printf("sum=%f jacobian_coeff=%f flux_param=%f\n", sum, jacobian_coeff, flux_param);

        for (int i = 0; i < helio_row_col.y; i++)
        {
            Vector<int> shadow_sectors;
            bool previous_true = false;
            for (int j = 0; j < helio_row_col.x; j++)
            {
                if (d_micro_sb[i * helio_row_col.x + j] && !previous_true)
                {
                    shadow_sectors.push_back(i);
                    previous_true = true;
                }
                if (!d_micro_sb[i * helio_row_col.x + j] && previous_true)
                {
                    shadow_sectors.push_back(i);
                    previous_true = false;
                }
            }
            if (d_micro_sb[i * helio_row_col.x + (helio_row_col.x - 1)])
            {
                shadow_sectors.push_back(helio_row_col.x - 1);
            }
            for (int j = 0; j < shadow_sectors.size(); j += 2)
            {
                float3 left_up = helio_vert0 + delta_x * shadow_sectors[j] / helio_row_col.x +
                                 delta_y * i / helio_row_col.y - focus_center;
                float3 right_up = helio_vert0 + delta_x * shadow_sectors[j + 1] / helio_row_col.x +
                                  delta_y * i / helio_row_col.y - focus_center;
                float3 right_bottom = helio_vert0 +
                                      delta_x * shadow_sectors[j + 1] / helio_row_col.x +
                                      delta_y * (i + 1) / helio_row_col.y - focus_center;
                float3 left_bottom = helio_vert0 + delta_x * shadow_sectors[j] / helio_row_col.x +
                                     delta_y * (i + 1) / helio_row_col.y - focus_center;
                // 左上 左下 右下 右上
                float2 new_micro_helio[4] = {
                    {dot(x_axis,      left_up), dot(y_axis,      left_up)},
                    {dot(x_axis,  left_bottom), dot(y_axis,  left_bottom)},
                    {dot(x_axis, right_bottom), dot(y_axis, right_bottom)},
                    {dot(x_axis,     right_up), dot(y_axis,     right_up)}
                };
                x13 = new_micro_helio[2].x;
                x10 = new_micro_helio[1].x;
                y11 = new_micro_helio[1].y;
                y10 = new_micro_helio[0].y;
                sum -= std::abs(1.0F * 0.25F *
                                (fast_math::erff((x13 - u) * sigma1_denominator) -
                                 fast_math::erff((x10 - u) * sigma1_denominator)) *
                                (fast_math::erff((y11 - v) * sigma2_denominator) -
                                 fast_math::erff((y10 - v) * sigma2_denominator)) *
                                jacobian_coeff * flux_param);
            }
        }
        // printf("%d\n", sum);

        sum *= project_cos;
        atomicAdd(&(d_array[my_id]), sum);
    }

    __global__ void convPlane(float* d_array, float sigma1, float sigma2, float flux_param,
                              float3 focus_center, int2 recv_row_col, float3 recv_vertex_0,
                              float3 d_x, float3 d_y, float3 plane_normal, int2 helio_row_col,
                              float3 helio_vert0, float3 delta_x, float3 delta_y, float3 x_axis,
                              float3 y_axis, float3 z_axis, float* gl_params, bool* d_micro_sb,
                              solar::ConvolutionModelType mode)
    {
        std::int64_t my_id = getThreadID();
        //    int2 row_col = hflcal->getRowCol();
        if (my_id >= recv_row_col.x * recv_row_col.y)
        {
            return;
        }
        // 计算当前点对应接收器平面中网格点坐标i,j
        int2 pos_ij = make_int2(my_id % recv_row_col.x, my_id / recv_row_col.x);
        // 计算网格点的世界坐标
        float3 pos = recv_vertex_0 + pos_ij.x * d_x + pos_ij.y * d_y;
        // 计算网格点在成像平面投影的世界坐标
        float3 projected_pos;
        // 斜平行投影找对应点
        std::tie(projected_pos, std::ignore) = intersect(z_axis, focus_center, pos, -z_axis);
        // 投影点在成像平面上的二维坐标
        float u = dot(x_axis, projected_pos - focus_center);
        float v = dot(y_axis, projected_pos - focus_center);

        float sum = 0;

        for (int i = 0; i < helio_row_col.x; i++)
        {
            for (int j = 0; j < helio_row_col.y; j++)
            {
                if (d_micro_sb[i * helio_row_col.x + j])
                {
                    continue;
                }
                float3 left_up = helio_vert0 + delta_x * i / helio_row_col.x +
                                 delta_y * j / helio_row_col.y - focus_center;
                float3 right_bottom = helio_vert0 + delta_x * (i + 1) / helio_row_col.x +
                                      delta_y * (j + 1) / helio_row_col.y - focus_center;
                float3 left_bottom = helio_vert0 + delta_x * i / helio_row_col.x +
                                     delta_y * (j + 1) / helio_row_col.y - focus_center;
                float xx = dot(left_up + right_bottom, x_axis) * 0.5F;
                float yy = dot(left_up + right_bottom, y_axis) * 0.5F;
                float3 tmp_x = left_bottom - left_up;
                float3 tmp_y = right_bottom - left_bottom;
                float2 delta_xx = make_float2(dot(tmp_x, x_axis), dot(tmp_x, y_axis)) * 0.5F;
                float2 delta_yy = make_float2(dot(tmp_y, x_axis), dot(tmp_y, y_axis)) * 0.5F;
                //            float delta_xx = dot(right_bottom-left_up, x_axis) * 0.5;
                //            float delta_yy = dot(right_bottom-left_up, y_axis) * 0.5;
                //            sum +=
                //            gauss_legendre(6,6,u,v,xx,yy,delta_xx,delta_yy,sigma1,sigma2,UNIZAR,flux_param,gl_params);
                //            sum +=
                //            gauss_legendre(6,6,u,v,xx,yy,delta_xx,delta_yy,sigma1,sigma2,Cauchy,flux_param,gl_params);
                sum += gaussLegendre(4, 4, u, v, xx, yy, delta_xx, delta_yy, sigma1, sigma2, mode,
                                     flux_param /*, gl_params*/);
            }
        }

        sum *= dot(z_axis, plane_normal);
        //    printf("%.3f\n",sum);
        //    if (dot(z_axis, plane_normal) >= 0)
        atomicAdd(&(d_array[my_id]), std::fmaxf(sum, 0.0F));
    }

    __global__ void convPlaneFast(float* d_array, float sigma1, float sigma2, float flux_param,
                                  float3 focus_center, int2 recv_row_col, float3 recv_vertex_0,
                                  float3 d_x, float3 d_y, float3 plane_normal, int2 helio_row_col,
                                  float3 x_axis, float3 y_axis, float3 z_axis, bool* d_micro_sb,
                                  float2* d_new_rectangle, float gradient_inv,
                                  solar::ConvolutionModelType mode)
    {
        std::int64_t my_id = getThreadID();
        //    int2 row_col = hflcal->getRowCol();
        if (my_id >= recv_row_col.x * recv_row_col.y)
        {
            return;
        }
        const auto project_cos = dot(z_axis, plane_normal);
        if (project_cos <= 0)
        {
            return;
        }
        // 计算当前点对应接收器平面中网格点坐标i,j
        const int2 pos_ij = make_int2(my_id % recv_row_col.x, my_id / recv_row_col.x);
        // 计算网格点的世界坐标
        const float3 pos = recv_vertex_0 + pos_ij.x * d_x + pos_ij.y * d_y;
        // 计算网格点在成像平面投影的世界坐标
        float3 projected_pos;
        // 斜平行投影找对应点
        std::tie(projected_pos, std::ignore) = intersect(z_axis, focus_center, pos, -z_axis);
        // 投影点在成像平面上的二维坐标
        const float u = dot(x_axis, projected_pos - focus_center);
        const float v = dot(y_axis, projected_pos - focus_center);
        // printf("%f, %f\n", u, v);

        const float x13 = d_new_rectangle[2].x;
        const float x10 = d_new_rectangle[1].x;
        const float y11 = d_new_rectangle[1].y;
        const float y10 = d_new_rectangle[0].y;
        const float sigma1_denominator = 1.0F / (std::numbers::sqrt2_v<float> * sigma1);
        const float sigma2_denominator = 1.0F / (std::numbers::sqrt2_v<float> * sigma2);
        const auto actual_u = u + gradient_inv * (v - y11);
        // const auto actual_u = u;
        float sum = std::abs(1.0F * 0.25F *
                             (fast_math::erff((x13 - actual_u) * sigma1_denominator) -
                              fast_math::erff((x10 - actual_u) * sigma1_denominator)) *
                             (fast_math::erff((y11 - v) * sigma2_denominator) -
                              fast_math::erff((y10 - v) * sigma2_denominator)) *
                             1.0F * flux_param);
        // printf("sum=%f jacobian_coeff=%f flux_param=%f\n", sum, jacobian_coeff, flux_param);
        // const auto aaa = fast_math::erff((x13 - actual_u) * sigma1_denominator);
        // const auto bbb = fast_math::erff((x10 - actual_u) * sigma1_denominator);
        // const auto ccc = fast_math::erff((y11 - v) * sigma2_denominator);
        // const auto ddd = fast_math::erff((y10 - v) * sigma2_denominator);
        // printf("%f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f\n", sigma1, sigma2, flux_param, x13, x10,
        // y11, y10, u, v, sum, aaa, bbb, ccc, ddd);

        const auto delta_x = std::abs(x13 - x10);
        const auto delta_y = std::abs(y10 - y11);

        for (int i = 0; i < helio_row_col.y; i++)
        {
            for (int j = 0; j < helio_row_col.x; j++)
            {
                if (!d_micro_sb[i * helio_row_col.x + j])
                {
                    continue;
                }
                const auto x13_1 = x10 + delta_x * ((j + 1) / helio_row_col.x);
                const auto x10_1 = x10 + delta_x * (j / helio_row_col.x);
                const auto y11_1 = y10 + delta_y * ((i + 1) / helio_row_col.y);
                const auto y10_1 = y10 + delta_y * (i / helio_row_col.y);
                const auto actual_u_1 = u + gradient_inv * (v - y11);
                sum -= std::abs(1.0F * 0.25F *
                                (fast_math::erff((x13_1 - actual_u_1) * sigma1_denominator) -
                                 fast_math::erff((x10_1 - actual_u_1) * sigma1_denominator)) *
                                (fast_math::erff((y11_1 - v) * sigma2_denominator) -
                                 fast_math::erff((y10_1 - v) * sigma2_denominator)) *
                                1.0F * flux_param);
            }
        }

        sum *= project_cos;
        //    printf("%.3f\n",sum);
        //    if (dot(z_axis, plane_normal) >= 0)
        atomicAdd(&(d_array[my_id]), sum);
    }

    __global__ void convPlaneFast2(float* d_array, dnn::UNIZAR2Output model_output,
                                   float3 focus_center, int2 recv_row_col, float3 recv_vertex_0,
                                   float3 d_x, float3 d_y, float3 plane_normal, int2 helio_row_col,
                                   float3 x_axis, float3 y_axis, float3 z_axis, bool* d_micro_sb,
                                   float gradient_inv, solar::ConvolutionModelType mode)
    {
        std::int64_t my_id = getThreadID();
        //    int2 row_col = hflcal->getRowCol();
        if (my_id >= recv_row_col.x * recv_row_col.y)
        {
            return;
        }
        const auto project_cos = dot(z_axis, plane_normal);
        if (project_cos <= 0)
        {
            return;
        }
        // 计算当前点对应接收器平面中网格点坐标i,j
        const int2 pos_ij = make_int2(my_id % recv_row_col.x, my_id / recv_row_col.x);
        // 计算网格点的世界坐标
        const float3 pos = recv_vertex_0 + pos_ij.x * d_x + pos_ij.y * d_y;
        // 计算网格点在成像平面投影的世界坐标
        float3 projected_pos;
        // 斜平行投影找对应点
        std::tie(projected_pos, std::ignore) = intersect(z_axis, focus_center, pos, -z_axis);
        // 投影点在成像平面上的二维坐标
        const float u = dot(x_axis, projected_pos - focus_center);
        const float v = dot(y_axis, projected_pos - focus_center);

        const float x13 = model_output.x3_;
        const float x10 = model_output.x0_;
        const float y11 = model_output.y1_;
        const float y10 = model_output.y0_;
        const float sigma1_denominator =
            1.0F / (std::numbers::sqrt2_v<float> * model_output.sigma_x_);
        const float sigma2_denominator =
            1.0F / (std::numbers::sqrt2_v<float> * model_output.sigma_y_);
        const auto actual_u = u + gradient_inv * (v - y11);
        float sum = std::abs(1.0F * 0.25F *
                             (fast_math::erff((x13 - actual_u) * sigma1_denominator) -
                              fast_math::erff((x10 - actual_u) * sigma1_denominator)) *
                             (fast_math::erff((y11 - v) * sigma2_denominator) -
                              fast_math::erff((y10 - v) * sigma2_denominator)) *
                             1.0F * model_output.peak_flux_param_);

        const auto delta_x = std::abs(x13 - x10);
        const auto delta_y = std::abs(y10 - y11);

        for (int i = 0; i < helio_row_col.y; i++)
        {
            for (int j = 0; j < helio_row_col.x; j++)
            {
                if (!d_micro_sb[i * helio_row_col.x + j])
                {
                    continue;
                }
                const auto x13_1 = x10 + delta_x * ((j + 1) / helio_row_col.x);
                const auto x10_1 = x10 + delta_x * (j / helio_row_col.x);
                const auto y11_1 = y10 + delta_y * ((i + 1) / helio_row_col.y);
                const auto y10_1 = y10 + delta_y * (i / helio_row_col.y);
                const auto actual_u_1 = u + gradient_inv * (v - y11);
                sum -= std::abs(1.0F * 0.25F *
                                (fast_math::erff((x13_1 - actual_u_1) * sigma1_denominator) -
                                 fast_math::erff((x10_1 - actual_u_1) * sigma1_denominator)) *
                                (fast_math::erff((y11_1 - v) * sigma2_denominator) -
                                 fast_math::erff((y10_1 - v) * sigma2_denominator)) *
                                1.0F * model_output.peak_flux_param_);
            }
        }

        sum *= project_cos;
        //    printf("%.3f\n",sum);
        //    if (dot(z_axis, plane_normal) >= 0)
        atomicAdd(&(d_array[my_id]), sum);
    }
} // namespace

ConvolutionMethod::ConvolutionMethod(float3* helio_verts, float3 helio_normal, float3 helio_pos,
                                     float3 recv_pos, float3 helio_size, float3 recv_size,
                                     float3 focus_center, int2 helio_row_col, int2 recv_row_col,
                                     bool* d_micro_sb, ConvolutionModelType mode, float* d_array,
                                     cudaStream_t stream)
    : d_micro_sb_(d_micro_sb), mode_(mode), recv_pos_(recv_pos),
      recv_size_(recv_size), helio_size_(helio_size), helio_normal_(helio_normal),
      focus_center_(focus_center), helio_pos_(helio_pos), recv_row_col_(recv_row_col),
      helio_row_col_(helio_row_col)
{
    if (stream != nullptr)
    {
        this->external_stream_ = true;
        this->stream_ = stream;
    }
    else
    {
        cudaStreamCreate(&stream_);
    }
    setHelioVerts(helio_verts);
    if (d_array != nullptr)
    {
        this->d_array_ = d_array;
        this->external_d_array_ = true;
    }
    else
    {
        CHECK_CUDA_ERRORS(cudaMallocAsync(
            &d_array_, sizeof(float) * recv_row_col.x * recv_row_col.y, this->stream_));
        CHECK_CUDA_ERRORS(cudaMemsetAsync(
            d_array_, 0, sizeof(float) * recv_row_col.x * recv_row_col.y, this->stream_));
    }

    // loadGlParam(10);
}

ConvolutionMethod::ConvolutionMethod(float3* helio_verts, float3 helio_normal,
                                     float3 helio_position, float3 focus_center, int2 helio_row_col,
                                     int2 recv_row_col, float3* recv_verts, float3 recv_normal,
                                     bool* d_micro_sb, float* d_array, cudaStream_t stream)
    : d_micro_sb_(d_micro_sb), helio_normal_(helio_normal), focus_center_(focus_center),
      helio_pos_(helio_position), plane_normal_(recv_normal), recv_row_col_(recv_row_col),
      helio_row_col_(helio_row_col)
{
    if (stream != nullptr)
    {
        this->external_stream_ = true;
        this->stream_ = stream;
    }
    else
    {
        cudaStreamCreate(&stream_);
    }
    setHelioVerts(helio_verts);
    std::copy(recv_verts, recv_verts + 4, recv_vertices_);
    if (d_array != nullptr)
    {
        this->d_array_ = d_array;
        this->external_d_array_ = true;
    }
    else
    {
        CHECK_CUDA_ERRORS(cudaMallocAsync(
            &d_array_, sizeof(float) * recv_row_col.x * recv_row_col.y, this->stream_));
        CHECK_CUDA_ERRORS(cudaMemsetAsync(
            d_array_, 0, sizeof(float) * recv_row_col.x * recv_row_col.y, this->stream_));
    }
}

void ConvolutionMethod::setHelioVerts(float3* const helio_vertices)
{
    std::copy(helio_vertices, helio_vertices + 4, helio_vertices_);
}

auto ConvolutionMethod::peekFitting(float peak, float3 x_axis, float3 y_axis, float3* interv,
                                    float area) -> float
{
    double left = 0.0F;
    double right = 6.0F;
    while (right - left > 0.001F)
    {
        double mid = (left + right) / 2;
        double val = intCauchyPeak(x_axis, y_axis, interv, mid, area);
        if (val > peak)
        {
            right = mid;
        }
        else
        {
            left = mid;
        }
    }
    return left;
}

auto ConvolutionMethod::intCauchyPeak(const float3& x_axis, const float3& y_axis,
                                      const float3* const interv, const float& s,
                                      const float& area) const -> float
{
    return solar::cpu::intCauchyPeak(structCopy<Float3>(x_axis), structCopy<Float3>(y_axis),
                                     std::bit_cast<const Float3* const>(interv), s, area,
                                     this->h_gl_params_.get(), structCopy<Float3>(focus_center_),
                                     sigma1_);
}

auto ConvolutionMethod::loadGlParam(int n) -> bool
{
    switch (n)
    {
    case 10:
        cpu2gpuAsync(d_gl_params_, gl_param_10_precomp, n * (n + 1), this->stream_);
        break;
    default:
        throw "Invalid gl-param num";
    }
    return true;
}

ConvolutionMethod::~ConvolutionMethod()
{
    if (d_array_ != nullptr && !this->external_d_array_)
    {
        CHECK_CUDA_ERRORS(cudaFreeAsync(d_array_, stream_));
        d_array_ = nullptr;
    }
    if (d_gl_params_ != nullptr)
    {
        CHECK_CUDA_ERRORS(cudaFreeAsync(d_gl_params_, stream_));
        d_gl_params_ = nullptr;
    }
    if (stream_ != nullptr && !this->external_stream_)
    {
        cudaStreamDestroy(stream_);
    }
}

void ConvolutionMethod::drawFluxMap()
{
    // 解析模型中定日镜镜面总能量计算
    // float distance = length(focus_center_ - helio_pos_);
    // float m_aa = (distance <= 1000)
    //                  ? (0.99321 - 0.0001176 * distance + 1.97 * 1e-8 * distance * distance)
    //                  : std::exp(-0.0001106 * distance);
    // // float s = helio_size_.x * helio_size_.z;
    // float cos_w = dot(normalize(focus_center_ - helio_pos_), helio_normal_);

    // 计算成像平面的x轴和y轴
    float3 x_axis;
    float3 y_axis;
    float3 z_axis;
    z_axis = normalize(helio_pos_ - focus_center_);
    float3 interv[3];
    for (int i = 0; i < 4; i++)
    {
        std::tie(interv[i], std::ignore) =
            intersect(z_axis, focus_center_, helio_vertices_[i], -z_axis);
    }
    x_axis = normalize(interv[1] - interv[0]);
    y_axis = cross(z_axis, x_axis);

    // float flux_param = 0;
    // if (flux_param < 1)
    // {
    //     flux_param = DNI_ * cos_w * helio_reflectivity_ * m_aa;
    // }

    // float area = length(cross(interv[2] - interv[1], interv[1] - interv[0]));
    // float p_w = length(interv[2] - interv[1]);
    // float p_h = area * cos_w / p_w;
    //    float tmp = 1;
    //    float tmp = sqrt(1-pow(z_axis.y,2));
    //    sigma1 = peekFitting_(peek / tmp, flux_param / area * 4.0*3.2, p_w, p_h);
    //    sigma2 = sigma1;
    //    flux_param *= fabs(dot(z_axis,plane_normal_));
    //    flux_param *=
    //    parallelProjectionParam(x_axis,y_axis,z_axis,focus_center_,helio_pos_,plane_normal_);

    // d_x, d_y
    float3 d_x = (recv_vertices_[1] - recv_vertices_[0]) / recv_row_col_.y;
    float3 d_y = (recv_vertices_[3] - recv_vertices_[0]) / recv_row_col_.x;
    float3 recv_vertex_0 = recv_vertices_[0] + d_x / 2 + d_y / 2;

    float3 delta_x = interv[2] - interv[1];
    float3 delta_y = interv[1] - interv[0];

    // 设置线程数
    int n_threads = 512;
    dim3 n_blocks;
    setThreadBlocks(n_blocks, n_threads, recv_row_col_.x * recv_row_col_.y, true);
    // debug check
    // float tmp = dot(z_axis, plane_normal_);
    int bias = helio_id_ * helio_row_col_.x * helio_row_col_.y;
    convPlane<<<n_blocks, n_threads, 0, stream_>>>(
        d_array_, sigma1_, sigma2_, flux_param_, focus_center_, recv_row_col_, recv_vertex_0, d_x,
        d_y, plane_normal_, helio_row_col_, interv[0], delta_x, delta_y, x_axis, y_axis, z_axis,
        d_gl_params_, &d_micro_sb_[bias], mode_);

    // cudaDeviceSynchronize();
    CHECK_CUDA_ERRORS(cudaGetLastError());
}

void ConvolutionMethod::drawFluxMapFast()
{
    // 计算成像平面的x轴和y轴
    float3 x_axis;
    float3 y_axis;
    float3 z_axis;
    z_axis = normalize(helio_pos_ - focus_center_);
    float3 interv[4];
    for (int i = 0; i < 4; i++)
    {
        std::tie(interv[i], std::ignore) =
            intersect(z_axis, focus_center_, helio_vertices_[i], -z_axis);
    }
    x_axis = normalize(interv[2] - interv[1]);
    y_axis = cross(z_axis, x_axis);

    // float3 delta_x = interv[2] - interv[1];
    // float3 delta_y = interv[1] - interv[0];

    float2 projected_interv[4];
    for (int i = 0; i < 4; i++)
    {
        projected_interv[i] = {dot(x_axis, interv[i] - focus_center_),
                               dot(y_axis, interv[i] - focus_center_)};
    }
    float2 current_projected_helio_verts[4];
    current_projected_helio_verts[0] = projected_interv[1];
    current_projected_helio_verts[1] = projected_interv[0];
    current_projected_helio_verts[2] = projected_interv[3];
    current_projected_helio_verts[3] = projected_interv[2];

    // const float x13 = current_projected_helio_verts[2].x;
    // const float x10 = current_projected_helio_verts[1].x;
    const float y11 = current_projected_helio_verts[1].y;
    const float y10 = current_projected_helio_verts[0].y;

    const auto projected_delta_y = y10 - y11;
    const auto projected_delta_x =
        current_projected_helio_verts[1].x - current_projected_helio_verts[0].x;
    const auto gradient_inv = projected_delta_x / projected_delta_y;

    // 左上 左下 右下 右上
    float2 new_rectangle[4] = {projected_interv[1], projected_interv[0], projected_interv[3],
                               projected_interv[2]};
    float2* d_new_rectangle = nullptr;
    cpu2gpuAsync(d_new_rectangle, new_rectangle, 4, stream_);

    // d_x, d_y
    float3 d_y = (recv_vertices_[1] - recv_vertices_[0]) / recv_row_col_.y;
    float3 d_x = (recv_vertices_[3] - recv_vertices_[0]) / recv_row_col_.x;
    float3 recv_vertex_0 = recv_vertices_[0] + d_x / 2 + d_y / 2;

    // float3 delta_x = interv[2] - interv[1];
    // float3 delta_y = interv[1] - interv[0];

    // 设置线程数
    int n_threads = 512;
    dim3 n_blocks;
    setThreadBlocks(n_blocks, n_threads, recv_row_col_.x * recv_row_col_.y, true);

    int bias = helio_id_ * helio_row_col_.x * helio_row_col_.y;
    if (!external_heliostat_projection_pos_)
    {
        convPlaneFast<<<n_blocks, n_threads, 0, stream_>>>(
            d_array_, sigma1_, sigma2_, flux_param_, focus_center_, recv_row_col_, recv_vertex_0,
            d_x, d_y, plane_normal_, helio_row_col_, x_axis, y_axis, z_axis, &d_micro_sb_[bias],
            d_new_rectangle, gradient_inv, mode_);
    }
    else
    {
        convPlaneFast2<<<n_blocks, n_threads, 0, stream_>>>(
            d_array_, {flux_param_, sigma1_, sigma2_, x3_, x0_, y1_, y0_}, focus_center_,
            recv_row_col_, recv_vertex_0, d_x, d_y, plane_normal_, helio_row_col_, x_axis, y_axis,
            z_axis, &d_micro_sb_[bias], gradient_inv, mode_);
    }

    // cudaDeviceSynchronize();
    CHECK_CUDA_ERRORS(cudaGetLastError());
}

void ConvolutionMethod::drawFluxMapCylinder()
{
    // 计算成像平面的x轴和y轴
    float3 x_axis;
    float3 y_axis;
    float3 z_axis;
    z_axis = normalize(helio_pos_ - focus_center_);
    float3 interv[4];
    for (int i = 0; i < 4; i++)
    {
        std::tie(interv[i], std::ignore) =
            intersect(z_axis, focus_center_, helio_vertices_[i], -z_axis);
    }
    x_axis = normalize(interv[2] - interv[1]);
    y_axis = cross(z_axis, x_axis);

    // 设置线程数
    int n_threads = 512;
    dim3 n_blocks;
    setThreadBlocks(n_blocks, n_threads, recv_row_col_.x * recv_row_col_.y, true);
    float3 cylinder_bottom_center =
        ::make_float3(recv_pos_.x, recv_pos_.y - 0.5f * recv_size_.y, recv_pos_.z);

    float3 delta_x = interv[2] - interv[1];
    float3 delta_y = interv[1] - interv[0];

    int bias = helio_id_ * helio_row_col_.x * helio_row_col_.y;

    convCylinder<<<n_blocks, n_threads, 0, stream_>>>(
        d_array_, 0, sigma1_, sigma2_, flux_param_, recv_row_col_, recv_size_,
        cylinder_bottom_center, helio_row_col_, interv[0], delta_x, delta_y, focus_center_, x_axis,
        y_axis, z_axis, d_gl_params_, &d_micro_sb_[bias], mode_);

    // cudaDeviceSynchronize();
    CHECK_CUDA_ERRORS(cudaGetLastError());
}

void ConvolutionMethod::drawFluxMapCylinderFast()
{
    // 计算成像平面的x轴和y轴
    float3 x_axis;
    float3 y_axis;
    float3 z_axis;
    z_axis = normalize(helio_pos_ - focus_center_);
    float3 interv[4];
    for (int i = 0; i < 4; i++)
    {
        std::tie(interv[i], std::ignore) =
            intersect(z_axis, focus_center_, helio_vertices_[i], -z_axis);
    }
    x_axis = normalize(interv[2] - interv[1]);
    y_axis = cross(z_axis, x_axis);

    // float3 delta_x = interv[2] - interv[1];
    // float3 delta_y = interv[1] - interv[0];

    float2 projected_interv[4];
    for (int i = 0; i < 4; i++)
    {
        projected_interv[i] = {dot(x_axis, interv[i] - focus_center_),
                               dot(y_axis, interv[i] - focus_center_)};
    }
    float2 current_projected_helio_verts[4];
    current_projected_helio_verts[0] = projected_interv[1];
    current_projected_helio_verts[1] = projected_interv[0];
    current_projected_helio_verts[2] = projected_interv[3];
    current_projected_helio_verts[3] = projected_interv[2];

    // const float x13 = current_projected_helio_verts[2].x;
    // const float x10 = current_projected_helio_verts[1].x;
    const float y11 = current_projected_helio_verts[1].y;
    const float y10 = current_projected_helio_verts[0].y;

    const auto projected_delta_y = y10 - y11;
    const auto projected_delta_x =
        current_projected_helio_verts[1].x - current_projected_helio_verts[0].x;
    const auto gradient_inv = projected_delta_x / projected_delta_y;

    // 左上 左下 右下 右上
    float2 new_rectangle[4] = {projected_interv[1], projected_interv[0], projected_interv[3],
                               projected_interv[2]};
    float jacobian_coeff = 1;
    float2* d_new_rectangle = nullptr;
    cpu2gpuAsync(d_new_rectangle, new_rectangle, 4, stream_);

    // 设置线程数
    int n_threads = 512;
    dim3 n_blocks;
    setThreadBlocks(n_blocks, n_threads, recv_row_col_.x * recv_row_col_.y, true);
    float3 cylinder_bottom_center =
        ::make_float3(recv_pos_.x, recv_pos_.y - 0.5f * recv_size_.y, recv_pos_.z);

    int bias = helio_id_ * helio_row_col_.x * helio_row_col_.y;

    if (!external_heliostat_projection_pos_)
    {
        convCylinderFast<<<n_blocks, n_threads, 0, stream_>>>(
            d_array_, 0, sigma1_, sigma2_, flux_param_, recv_row_col_, recv_size_,
            cylinder_bottom_center, helio_row_col_, focus_center_, x_axis, y_axis, z_axis,
            &d_micro_sb_[bias], d_new_rectangle, jacobian_coeff, gradient_inv, mode_);
    }
    else
    {
        convCylinderFast2<<<n_blocks, n_threads, 0, stream_>>>(
            d_array_, 0, {flux_param_, sigma1_, sigma2_, x3_, x0_, y1_, y0_}, recv_row_col_,
            recv_size_, cylinder_bottom_center, helio_row_col_, focus_center_, x_axis, y_axis,
            z_axis, &d_micro_sb_[bias], jacobian_coeff, gradient_inv, mode_);
    }

    // cudaDeviceSynchronize();
    CHECK_CUDA_ERRORS(cudaGetLastError());
}

void ConvolutionMethod::setCudaStream(const cudaStream_t& stream)
{
    // if (this->stream_ != nullptr && !this->external_stream_)
    // {
    //     cudaStreamDestroy(this->stream_);
    // }
    this->stream_ = stream;
}

struct alignas(16) vec3d
{
        float x, y, z;
};

ConvolutionMethod::ConvolutionMethod(float3* d_image_plane_xyz_axes, float3* d_helio_verts,
                                     bool* d_micro_sb, bool* d_micro_sb_per_row,
                                     float3* d_flux_params, float3* d_focus_centers,
                                     std::size_t num_helios, float3 recv_size,
                                     float3 recv_pos_center, int2 helio_row_col, int2 recv_row_col,
                                     ConvolutionModelType mode, float* d_array, cudaStream_t stream)
{
    if (mode == ConvolutionModelType::kiCauchy)
    {
        const auto error_str =
            solar::format("[ERROR] In function {}: Model iCauchy is not supported for now\n",
                          PRETTY_FUNCTION_NAME);
        std::cout << error_str;
        throw std::invalid_argument(error_str);
    }
    cudaDeviceSetCacheConfig(cudaFuncCachePreferL1);
    this->external_d_array_ = true;
    if (stream != nullptr)
    {
        this->external_stream_ = true;
        this->stream_ = stream;
    }
    else
    {
        CHECK_CUDA_ERRORS(cudaStreamCreate(&stream_));
    }
    // 接收器底部中心坐标
    const float3 recv_pos =
        ::make_float3(recv_pos_center.x, recv_pos_center.y - 0.5F * recv_size.y, recv_pos_center.z);
    const std::size_t num_helio_vertices = num_helios * 4;
    // 剪切变换后的定日镜顶点（成像平面上）
    float2* d_projected_helio_verts = nullptr;
    float* d_gradient_inv_per_heliostat = nullptr;
    CHECK_CUDA_ERRORS(cudaMallocAsync(std::addressof(d_projected_helio_verts),
                                      sizeof(float2) * num_helio_vertices, stream_));
    CHECK_CUDA_ERRORS(cudaMallocAsync(std::addressof(d_gradient_inv_per_heliostat),
                                      sizeof(float) * num_helios, stream_));
    // const auto focus_center = this->focus_center_;

    // 计算定日镜在成像平面上的投影的长宽、剪切变换后的定日镜顶点（成像平面上）
    thrust::async::for_each(
        thrust::device.on(stream_), thrust::make_counting_iterator(std::size_t(0)),
        thrust::make_counting_iterator(num_helios),
        [d_image_plane_xyz_axes, d_helio_verts, d_focus_centers,
         d_gradient_inv_per_heliostat, d_projected_helio_verts] __device__(const std::size_t& index)
        {
            const auto focus_center = d_focus_centers[index];
            const auto current_image_plane_axis_ptr = d_image_plane_xyz_axes + (index * 3);
            const auto current_image_plane_x_axis = current_image_plane_axis_ptr[0];
            const auto current_image_plane_y_axis = current_image_plane_axis_ptr[1];
            // const auto current_image_plane_z_axis = current_image_plane_axis_ptr[2];
            const auto current_helio_verts_ptr = d_helio_verts + (index * 4);
            // float3 interv[4];
            // for (int i = 0; i < 4; i++)
            // {
            //     std::tie(interv[i], std::ignore) =
            //         intersect(current_image_plane_z_axis, focus_center, current_helio_verts_ptr[i],
            //                   -current_image_plane_z_axis);
            // }

            float2 projected_interv[4];
            for (int i = 0; i < 4; i++)
            {
                projected_interv[i] = {dot(current_image_plane_x_axis, current_helio_verts_ptr[i] - focus_center),
                                       dot(current_image_plane_y_axis, current_helio_verts_ptr[i] - focus_center)};
            }
            // 左上 左下 右下 右上
            auto current_projected_helio_verts = d_projected_helio_verts + (index * 4);
            current_projected_helio_verts[0] = projected_interv[1];
            current_projected_helio_verts[1] = projected_interv[0];
            current_projected_helio_verts[2] = projected_interv[3];
            current_projected_helio_verts[3] = projected_interv[2];

            // const float x13 = current_projected_helio_verts[2].x;
            // const float x10 = current_projected_helio_verts[1].x;
            const float y11 = current_projected_helio_verts[1].y;
            const float y10 = current_projected_helio_verts[0].y;

            const auto projected_delta_y = y10 - y11;
            const auto projected_delta_x =
                current_projected_helio_verts[1].x - current_projected_helio_verts[0].x;
            const auto gradient_inv = projected_delta_x / projected_delta_y;
            d_gradient_inv_per_heliostat[index] = gradient_inv;
            // printf("x13=%f x10=%f y11=%f y10=%f\n", x13, x10, y11, y10);
        });

    const auto recv_row_col_mul = recv_row_col.x * recv_row_col.y;
    const float theta_0 = 0.0F;

    float3* d_receiver_pixel_pos = nullptr;
    float* d_pixel_angles = nullptr;
    CHECK_CUDA_ERRORS(
        cudaMallocAsync(std::addressof(d_pixel_angles), sizeof(float) * recv_row_col.x, stream_));
    CHECK_CUDA_ERRORS(cudaMallocAsync(std::addressof(d_receiver_pixel_pos),
                                      sizeof(float3) * recv_row_col_mul, stream_));
    const float delta_pixel_angle = 1.0F / recv_row_col.x * 2.0F * inter_op::kMathPi;
    float pixel_angles_0[1] = {theta_0 + 0.5F * delta_pixel_angle};
    cpu2gpuAsync(d_pixel_angles, pixel_angles_0, 1);
    thrust::async::for_each(
        thrust::device.on(stream_), thrust::make_counting_iterator(std::size_t(1)),
        thrust::make_counting_iterator(std::size_t(recv_row_col.x)),
        [d_pixel_angles, delta_pixel_angle] __device__(const std::size_t& index)
        { d_pixel_angles[index] = d_pixel_angles[0] + delta_pixel_angle * index; });

    // 计算接收器上每个像素的世界坐标
    thrust::async::for_each(
        thrust::device.on(stream_), thrust::make_counting_iterator(std::size_t(0)),
        thrust::make_counting_iterator(std::size_t(recv_row_col_mul)),
        [d_receiver_pixel_pos, recv_pos, recv_size, recv_row_col,
         d_pixel_angles] __device__(const std::size_t& my_id)
        {
            // 计算当前点对应接收器平面中网格点坐标i,j
            const int2 pos_ij = ::make_int2(my_id % recv_row_col.x, my_id / recv_row_col.x);
            // 计算网格点的世界坐标
            const float pixel_angle = d_pixel_angles[pos_ij.x];
            d_receiver_pixel_pos[my_id] =
                recv_pos +
                ::make_float3(__cosf(pixel_angle) * recv_size.x, 0,
                              __sinf(pixel_angle) * recv_size.x) +
                ::make_float3(0, 1.0F * (pos_ij.y + 0.5F) / recv_row_col.y * recv_size.y, 0);
        });

    const std::size_t helio_row_col_mul = helio_row_col.x * helio_row_col.y;
    // constexpr std::uint16_t bool_true_2 = 0x101;
    // constexpr std::uint32_t bool_true_4 = 0x1010101;

    // 整个镜场仿真（卷积模型 exp-based）
    // 微定日镜划分较粗，不合并计算积分
    thrust::async::for_each(
        thrust::device.on(stream_), thrust::make_counting_iterator(std::size_t(0)),
        thrust::make_counting_iterator(std::size_t(num_helios * recv_row_col_mul)),
        [d_image_plane_xyz_axes, d_flux_params, recv_row_col,
         recv_row_col_mul, d_focus_centers, helio_row_col, d_micro_sb, d_micro_sb_per_row, d_array,
         d_receiver_pixel_pos, d_pixel_angles, d_gradient_inv_per_heliostat, d_projected_helio_verts,
         helio_row_col_mul] __device__(const std::size_t& index)
        {
            const auto helio_id = index / (recv_row_col_mul);
            const auto my_id = index % (recv_row_col_mul);
            // printf("helio_id=%d\n", helio_id);
            // // 计算当前点对应接收器平面中网格点坐标i,j
            const int2 pos_ij = ::make_int2(my_id % recv_row_col.x, my_id / recv_row_col.x);
            // // 计算网格点的世界坐标
            // const float pixel_angle =
            //     theta_0 + 1.0F * (pos_ij.x + 0.5F) / recv_row_col.x * 2.0F *
            //     inter_op::kMathPi;
            auto* const current_image_plane_axis_ptr = d_image_plane_xyz_axes + (helio_id * 3);
            const auto x_axis = current_image_plane_axis_ptr[0];
            const auto y_axis = current_image_plane_axis_ptr[1];
            const auto z_axis = current_image_plane_axis_ptr[2];
            // printf("%f %f %f\n", x_axis.x, x_axis.y, x_axis.z);
            const float pixel_angle = d_pixel_angles[pos_ij.x];
            const auto project_cos =
                dot(z_axis, ::make_float3(__cosf(pixel_angle), 0, __sinf(pixel_angle)));
            if (project_cos <= 0)
            {
                return;
            }
            const float3 pos = d_receiver_pixel_pos[my_id];
            // 计算网格点在成像平面投影的世界坐标
            // float3 projected_pos;
            // 斜平行投影找对应点
            const auto focus_center = d_focus_centers[helio_id];
            // std::tie(projected_pos, std::ignore) = intersect(z_axis, focus_center, pos, -z_axis);
            // 投影点在成像平面上的二维坐标
            const float u = dot(x_axis, pos - focus_center); // x0
            const float v = dot(y_axis, pos - focus_center); // y0

            const auto current_flux_params = d_flux_params[helio_id];
            const auto sigma1 = current_flux_params.x;
            const auto sigma2 = current_flux_params.y;
            const auto flux_param = current_flux_params.z;

            auto* const current_projected_helio_verts = d_projected_helio_verts + (helio_id * 4);
            const float x13 = current_projected_helio_verts[2].x;
            const float x10 = current_projected_helio_verts[1].x;
            const float y11 = current_projected_helio_verts[1].y;
            const float y10 = current_projected_helio_verts[0].y;
            const auto current_gradient_inv = d_gradient_inv_per_heliostat[helio_id];
            const auto actual_u = u + current_gradient_inv * (v - y11);
            const float sigma1_denominator = 1.0F / (std::numbers::sqrt2_v<float> * sigma1);
            const float sigma2_denominator = 1.0F / (std::numbers::sqrt2_v<float> * sigma2);
            float sum = std::abs(1.0F * 0.25F *
                                 (fast_math::erff((x13 - actual_u) * sigma1_denominator) -
                                  fast_math::erff((x10 - actual_u) * sigma1_denominator)) *
                                 (fast_math::erff((y11 - v) * sigma2_denominator) -
                                  fast_math::erff((y10 - v) * sigma2_denominator)) *
                                 1.0F /*jacobian_coeff*/ * flux_param);
            // printf("xxx %f, %f, %f, %f, %f, %f, %f, %f, %f, %f\n", sigma1, sigma2,
            // flux_param, x13, x10, y11, y10, u, v, sum);

            const auto delta_x = std::abs(x13 - x10);
            const auto delta_y = std::abs(y10 - y11);
            auto* const current_micro_sb = d_micro_sb + helio_id * helio_row_col_mul;
            auto* const current_micro_sb_per_row = d_micro_sb_per_row + helio_id * helio_row_col.y;
            using wide_read_type = std::uint32_t;
            constexpr std::size_t mask = ~(sizeof(wide_read_type) - 1);
            constexpr auto for_step = sizeof(wide_read_type) / sizeof(bool);

            for (int row_index1 = 0; row_index1 < helio_row_col.y; row_index1 += for_step)
            {
                const auto* current_micro_sb_per_row_ptr =
                    std::addressof(current_micro_sb_per_row[row_index1]);
                wide_read_type current_micro_sb_per_row_wide(0);
                std::memcpy(std::addressof(current_micro_sb_per_row_wide),
                            current_micro_sb_per_row_ptr, sizeof(wide_read_type));
                if (current_micro_sb_per_row_wide != 0)
                {
                    for (std::size_t row_index2 = 0; row_index2 < for_step; row_index2++)
                    {
                        const auto current_row_index = row_index1 + row_index2;
                        if (current_micro_sb_per_row[current_row_index] != 0)
                        {
                            auto* const current_micro_sb_row = std::addressof(
                                current_micro_sb[current_row_index * helio_row_col.x]);
                            for (std::size_t index1 = 0; index1 < (helio_row_col.x & mask);
                                 index1 += for_step)
                            {
                                const auto* current_micro_sb_ptr =
                                    std::addressof(current_micro_sb_row[index1]);
                                wide_read_type current_micro_sb_wide(0);
                                std::memcpy(std::addressof(current_micro_sb_wide),
                                            current_micro_sb_ptr, sizeof(wide_read_type));
                                if (current_micro_sb_wide != 0)
                                {
                                    // #pragma unroll
                                    for (std::size_t index2 = 0; index2 < for_step; index2++)
                                    {
                                        // const auto i = (index1 + index2) / helio_row_col.x;
                                        const auto j = index1 + index2;
                                        if (!current_micro_sb[row_index1 * helio_row_col.x + j])
                                        {
                                            continue;
                                        }
                                        const auto x13_1 =
                                            x10 + delta_x * ((j + 1) / helio_row_col.x);
                                        const auto x10_1 = x10 + delta_x * (j / helio_row_col.x);
                                        const auto y11_1 =
                                            y10 + delta_y * ((row_index1 + 1) / helio_row_col.y);
                                        const auto y10_1 =
                                            y10 + delta_y * (row_index1 / helio_row_col.y);
                                        const auto actual_u_1 =
                                            u + current_gradient_inv * (v - y11);
                                        sum -= std::abs(
                                            1.0F * 0.25F *
                                            (fast_math::erff((x13_1 - actual_u_1) *
                                                             sigma1_denominator) -
                                             fast_math::erff((x10_1 - actual_u_1) *
                                                             sigma1_denominator)) *
                                            (fast_math::erff((y11_1 - v) * sigma2_denominator) -
                                             fast_math::erff((y10_1 - v) * sigma2_denominator)) *
                                            1.0F /*jacobian_coeff*/ * flux_param);
                                    }
                                    for (std::size_t index1 = (helio_row_col.x & mask);
                                         index1 < helio_row_col.x; index1++)
                                    {
                                        // const auto i = index1 / helio_row_col.x;
                                        const auto j = index1;
                                        if (!current_micro_sb[row_index1 * helio_row_col.x + j])
                                        {
                                            continue;
                                        }
                                        const auto x13_1 =
                                            x10 + delta_x * ((j + 1) / helio_row_col.x);
                                        const auto x10_1 = x10 + delta_x * (j / helio_row_col.x);
                                        const auto y11_1 =
                                            y10 + delta_y * ((row_index1 + 1) / helio_row_col.y);
                                        const auto y10_1 =
                                            y10 + delta_y * (row_index1 / helio_row_col.y);
                                        const auto actual_u_1 =
                                            u + current_gradient_inv * (v - y11);
                                        sum -= std::abs(
                                            1.0F * 0.25F *
                                            (fast_math::erff((x13_1 - actual_u_1) *
                                                             sigma1_denominator) -
                                             fast_math::erff((x10_1 - actual_u_1) *
                                                             sigma1_denominator)) *
                                            (fast_math::erff((y11_1 - v) * sigma2_denominator) -
                                             fast_math::erff((y10_1 - v) * sigma2_denominator)) *
                                            1.0F /*jacobian_coeff*/ * flux_param);
                                    }
                                }
                            }
                        }
                    }
                }
            }

            sum *= project_cos;
            // printf("%f, %f, %f, %f\n", sigma1, sigma2, flux_param, sum);
            atomicAdd(std::addressof(d_array[my_id]), sum);
        });

    CHECK_CUDA_ERRORS(cudaFreeAsync(d_projected_helio_verts, stream_));
    CHECK_CUDA_ERRORS(cudaFreeAsync(d_receiver_pixel_pos, stream_));
    CHECK_CUDA_ERRORS(cudaFreeAsync(d_pixel_angles, stream_));
    if (!external_stream_)
    {
        CHECK_CUDA_ERRORS(cudaStreamSynchronize(stream_));
    }
}

ConvolutionMethod::ConvolutionMethod(float3* d_image_plane_xyz_axes, float3* d_helio_verts,
                                     bool* d_micro_sb, bool* d_micro_sb_per_row,
                                     ::solar::dnn::UNIZAR2Output* d_model_outputs,
                                     float3* d_focus_centers, std::size_t num_helios,
                                     float3 recv_size, float3 recv_pos_center, int2 helio_row_col,
                                     int2 recv_row_col, ConvolutionModelType mode, float* d_array,
                                     cudaStream_t stream)
{
    if (mode == ConvolutionModelType::kiCauchy)
    {
        const auto error_str =
            solar::format("[ERROR] In function {}: Model iCauchy is not supported for now\n",
                          PRETTY_FUNCTION_NAME);
        std::cout << error_str;
        throw std::invalid_argument(error_str);
    }
    cudaDeviceSetCacheConfig(cudaFuncCachePreferL1);
    this->external_d_array_ = true;
    if (stream != nullptr)
    {
        this->external_stream_ = true;
        this->stream_ = stream;
    }
    else
    {
        CHECK_CUDA_ERRORS(cudaStreamCreate(&stream_));
    }
    // 接收器底部中心坐标
    const float3 recv_pos =
        ::make_float3(recv_pos_center.x, recv_pos_center.y - 0.5F * recv_size.y, recv_pos_center.z);
    // const std::size_t num_helio_vertices = num_helios * 4;
    // 剪切变换后的定日镜顶点（成像平面上）
    // float2* d_projected_helio_verts = nullptr;
    float* d_gradient_inv_per_heliostat = nullptr;
    // CHECK_CUDA_ERRORS(cudaMallocAsync(std::addressof(d_projected_helio_verts),
    //                                   sizeof(float2) * num_helio_vertices, stream_));
    CHECK_CUDA_ERRORS(cudaMallocAsync(std::addressof(d_gradient_inv_per_heliostat),
                                      sizeof(float) * num_helios, stream_));
    // const auto focus_center = this->focus_center_;

    // 计算定日镜在成像平面上的投影的长宽、剪切变换后的定日镜顶点（成像平面上）
    thrust::async::for_each(
        thrust::device.on(stream_), thrust::make_counting_iterator(std::size_t(0)),
        thrust::make_counting_iterator(num_helios),
        [d_image_plane_xyz_axes, d_helio_verts, d_focus_centers,
         d_gradient_inv_per_heliostat] __device__(const std::size_t& index)
        {
            const auto focus_center = d_focus_centers[index];
            const auto current_image_plane_axis_ptr = d_image_plane_xyz_axes + (index * 3);
            const auto current_image_plane_x_axis = current_image_plane_axis_ptr[0];
            const auto current_image_plane_y_axis = current_image_plane_axis_ptr[1];
            // const auto current_image_plane_z_axis = current_image_plane_axis_ptr[2];
            const auto current_helio_verts_ptr = d_helio_verts + (index * 4);
            // float3 interv[4];
            // for (int i = 0; i < 4; i++)
            // {
            //     std::tie(interv[i], std::ignore) =
            //         intersect(current_image_plane_z_axis, focus_center, current_helio_verts_ptr[i],
            //                   -current_image_plane_z_axis);
            // }

            float2 projected_interv[4];
            for (int i = 0; i < 4; i++)
            {
                projected_interv[i] = {dot(current_image_plane_x_axis, current_helio_verts_ptr[i] - focus_center),
                                       dot(current_image_plane_y_axis, current_helio_verts_ptr[i] - focus_center)};
            }
            // 左上 左下 右下 右上
            float2 current_projected_helio_verts[2];
            current_projected_helio_verts[0] = projected_interv[1];
            current_projected_helio_verts[1] = projected_interv[0];
            // current_projected_helio_verts[2] = projected_interv[3];
            // current_projected_helio_verts[3] = projected_interv[2];

            // const float x13 = current_projected_helio_verts[2].x;
            // const float x10 = current_projected_helio_verts[1].x;
            const float y11 = current_projected_helio_verts[1].y;
            const float y10 = current_projected_helio_verts[0].y;

            const auto projected_delta_y = y10 - y11;
            const auto projected_delta_x =
                current_projected_helio_verts[1].x - current_projected_helio_verts[0].x;
            const auto gradient_inv = projected_delta_x / projected_delta_y;
            d_gradient_inv_per_heliostat[index] = gradient_inv;
            // printf("x13=%f x10=%f y11=%f y10=%f\n", x13, x10, y11, y10);
        });

    const auto recv_row_col_mul = recv_row_col.x * recv_row_col.y;
    const float theta_0 = 0.0F;

    float3* d_receiver_pixel_pos = nullptr;
    float* d_pixel_angles = nullptr;
    CHECK_CUDA_ERRORS(
        cudaMallocAsync(std::addressof(d_pixel_angles), sizeof(float) * recv_row_col.x, stream_));
    CHECK_CUDA_ERRORS(cudaMallocAsync(std::addressof(d_receiver_pixel_pos),
                                      sizeof(float3) * recv_row_col_mul, stream_));
    const float delta_pixel_angle = 1.0F / recv_row_col.x * 2.0F * inter_op::kMathPi;
    float pixel_angles_0[1] = {theta_0 + 0.5F * delta_pixel_angle};
    cpu2gpuAsync(d_pixel_angles, pixel_angles_0, 1);
    thrust::async::for_each(
        thrust::device.on(stream_), thrust::make_counting_iterator(std::size_t(1)),
        thrust::make_counting_iterator(std::size_t(recv_row_col.x)),
        [d_pixel_angles, delta_pixel_angle] __device__(const std::size_t& index)
        { d_pixel_angles[index] = d_pixel_angles[0] + delta_pixel_angle * index; });

    // 计算接收器上每个像素的世界坐标
    thrust::async::for_each(
        thrust::device.on(stream_), thrust::make_counting_iterator(std::size_t(0)),
        thrust::make_counting_iterator(std::size_t(recv_row_col_mul)),
        [d_receiver_pixel_pos, recv_pos, recv_size, recv_row_col,
         d_pixel_angles] __device__(const std::size_t& my_id)
        {
            // 计算当前点对应接收器平面中网格点坐标i,j
            const int2 pos_ij = ::make_int2(my_id % recv_row_col.x, my_id / recv_row_col.x);
            // 计算网格点的世界坐标
            const float pixel_angle = d_pixel_angles[pos_ij.x];
            d_receiver_pixel_pos[my_id] =
                recv_pos +
                ::make_float3(__cosf(pixel_angle) * recv_size.x, 0,
                              __sinf(pixel_angle) * recv_size.x) +
                ::make_float3(0, 1.0F * (pos_ij.y + 0.5F) / recv_row_col.y * recv_size.y, 0);
        });

    const std::size_t helio_row_col_mul = helio_row_col.x * helio_row_col.y;
    // constexpr std::uint16_t bool_true_2 = 0x101;
    // constexpr std::uint32_t bool_true_4 = 0x1010101;

    // 整个镜场仿真（卷积模型 exp-based）
    // 微定日镜划分较粗，不合并计算积分
    thrust::async::for_each(
        thrust::device.on(stream_), thrust::make_counting_iterator(std::size_t(0)),
        thrust::make_counting_iterator(std::size_t(num_helios * recv_row_col_mul)),
        [d_image_plane_xyz_axes, d_model_outputs, recv_row_col, recv_row_col_mul, d_focus_centers,
         helio_row_col, d_micro_sb, d_micro_sb_per_row, d_array, d_receiver_pixel_pos,
         d_pixel_angles, d_gradient_inv_per_heliostat,
         helio_row_col_mul] __device__(const std::size_t& index)
        {
            const auto helio_id = index / (recv_row_col_mul);
            const auto my_id = index % (recv_row_col_mul);
            // printf("helio_id=%d\n", helio_id);
            // // 计算当前点对应接收器平面中网格点坐标i,j
            const int2 pos_ij = ::make_int2(my_id % recv_row_col.x, my_id / recv_row_col.x);
            // // 计算网格点的世界坐标
            // const float pixel_angle =
            //     theta_0 + 1.0F * (pos_ij.x + 0.5F) / recv_row_col.x * 2.0F *
            //     inter_op::kMathPi;
            auto* const current_image_plane_axis_ptr = d_image_plane_xyz_axes + (helio_id * 3);
            const auto x_axis = current_image_plane_axis_ptr[0];
            const auto y_axis = current_image_plane_axis_ptr[1];
            const auto z_axis = current_image_plane_axis_ptr[2];
            // printf("%f %f %f\n", x_axis.x, x_axis.y, x_axis.z);
            const float pixel_angle = d_pixel_angles[pos_ij.x];
            const auto project_cos =
                dot(z_axis, ::make_float3(__cosf(pixel_angle), 0, __sinf(pixel_angle)));
            if (project_cos <= 0)
            {
                return;
            }
            const float3 pos = d_receiver_pixel_pos[my_id];
            // 计算网格点在成像平面投影的世界坐标
            // float3 projected_pos;
            // 斜平行投影找对应点
            const auto focus_center = d_focus_centers[helio_id];
            // std::tie(projected_pos, std::ignore) = intersect(z_axis, focus_center, pos, -z_axis);
            // 投影点在成像平面上的二维坐标
            const float u = dot(x_axis, pos - focus_center); // x0
            const float v = dot(y_axis, pos - focus_center); // y0

            const auto current_model_params = d_model_outputs[helio_id];
            const auto sigma1 = current_model_params.sigma_x_;
            const auto sigma2 = current_model_params.sigma_y_;
            const auto flux_param = current_model_params.peak_flux_param_;

            const float x13 = current_model_params.x3_;
            const float x10 = current_model_params.x0_;
            const float y11 = current_model_params.y1_;
            const float y10 = current_model_params.y0_;
            const auto current_gradient_inv = d_gradient_inv_per_heliostat[helio_id];
            const auto actual_u = u + current_gradient_inv * (v - y11);
            const float sigma1_denominator = 1.0F / (std::numbers::sqrt2_v<float> * sigma1);
            const float sigma2_denominator = 1.0F / (std::numbers::sqrt2_v<float> * sigma2);
            float sum = std::abs(1.0F * 0.25F *
                                 (fast_math::erff((x13 - actual_u) * sigma1_denominator) -
                                  fast_math::erff((x10 - actual_u) * sigma1_denominator)) *
                                 (fast_math::erff((y11 - v) * sigma2_denominator) -
                                  fast_math::erff((y10 - v) * sigma2_denominator)) *
                                 1.0F /*jacobian_coeff*/ * flux_param);
            // printf("xxx %f, %f, %f, %f, %f, %f, %f, %f, %f, %f\n", sigma1, sigma2,
            // flux_param, x13, x10, y11, y10, u, v, sum);

            const auto delta_x = std::abs(x13 - x10);
            const auto delta_y = std::abs(y10 - y11);
            auto* const current_micro_sb = d_micro_sb + helio_id * helio_row_col_mul;
            auto* const current_micro_sb_per_row = d_micro_sb_per_row + helio_id * helio_row_col.y;
            using wide_read_type = std::uint32_t;
            constexpr std::size_t mask = ~(sizeof(wide_read_type) - 1);
            constexpr auto for_step = sizeof(wide_read_type) / sizeof(bool);

            for (int row_index1 = 0; row_index1 < helio_row_col.y; row_index1 += for_step)
            {
                const auto* current_micro_sb_per_row_ptr =
                    std::addressof(current_micro_sb_per_row[row_index1]);
                wide_read_type current_micro_sb_per_row_wide(0);
                std::memcpy(std::addressof(current_micro_sb_per_row_wide),
                            current_micro_sb_per_row_ptr, sizeof(wide_read_type));
                if (current_micro_sb_per_row_wide != 0)
                {
                    for (std::size_t row_index2 = 0; row_index2 < for_step; row_index2++)
                    {
                        const auto current_row_index = row_index1 + row_index2;
                        if (current_micro_sb_per_row[current_row_index] != 0)
                        {
                            auto* const current_micro_sb_row = std::addressof(
                                current_micro_sb[current_row_index * helio_row_col.x]);
                            for (std::size_t index1 = 0; index1 < (helio_row_col.x & mask);
                                 index1 += for_step)
                            {
                                const auto* current_micro_sb_ptr =
                                    std::addressof(current_micro_sb_row[index1]);
                                wide_read_type current_micro_sb_wide(0);
                                std::memcpy(std::addressof(current_micro_sb_wide),
                                            current_micro_sb_ptr, sizeof(wide_read_type));
                                if (current_micro_sb_wide != 0)
                                {
                                    // #pragma unroll
                                    for (std::size_t index2 = 0; index2 < for_step; index2++)
                                    {
                                        // const auto i = (index1 + index2) / helio_row_col.x;
                                        const auto j = index1 + index2;
                                        if (!current_micro_sb[row_index1 * helio_row_col.x + j])
                                        {
                                            continue;
                                        }
                                        const auto x13_1 =
                                            x10 + delta_x * ((j + 1) / helio_row_col.x);
                                        const auto x10_1 = x10 + delta_x * (j / helio_row_col.x);
                                        const auto y11_1 =
                                            y10 + delta_y * ((row_index1 + 1) / helio_row_col.y);
                                        const auto y10_1 =
                                            y10 + delta_y * (row_index1 / helio_row_col.y);
                                        const auto actual_u_1 =
                                            u + current_gradient_inv * (v - y11);
                                        sum -= std::abs(
                                            1.0F * 0.25F *
                                            (fast_math::erff((x13_1 - actual_u_1) *
                                                             sigma1_denominator) -
                                             fast_math::erff((x10_1 - actual_u_1) *
                                                             sigma1_denominator)) *
                                            (fast_math::erff((y11_1 - v) * sigma2_denominator) -
                                             fast_math::erff((y10_1 - v) * sigma2_denominator)) *
                                            1.0F /*jacobian_coeff*/ * flux_param);
                                    }
                                    for (std::size_t index1 = (helio_row_col.x & mask);
                                         index1 < helio_row_col.x; index1++)
                                    {
                                        // const auto i = index1 / helio_row_col.x;
                                        const auto j = index1;
                                        if (!current_micro_sb[row_index1 * helio_row_col.x + j])
                                        {
                                            continue;
                                        }
                                        const auto x13_1 =
                                            x10 + delta_x * ((j + 1) / helio_row_col.x);
                                        const auto x10_1 = x10 + delta_x * (j / helio_row_col.x);
                                        const auto y11_1 =
                                            y10 + delta_y * ((row_index1 + 1) / helio_row_col.y);
                                        const auto y10_1 =
                                            y10 + delta_y * (row_index1 / helio_row_col.y);
                                        const auto actual_u_1 =
                                            u + current_gradient_inv * (v - y11);
                                        sum -= std::abs(
                                            1.0F * 0.25F *
                                            (fast_math::erff((x13_1 - actual_u_1) *
                                                             sigma1_denominator) -
                                             fast_math::erff((x10_1 - actual_u_1) *
                                                             sigma1_denominator)) *
                                            (fast_math::erff((y11_1 - v) * sigma2_denominator) -
                                             fast_math::erff((y10_1 - v) * sigma2_denominator)) *
                                            1.0F /*jacobian_coeff*/ * flux_param);
                                    }
                                }
                            }
                        }
                    }
                }
            }

            sum *= project_cos;
            // printf("%f, %f, %f, %f\n", sigma1, sigma2, flux_param, sum);
            atomicAdd(std::addressof(d_array[my_id]), sum);
        });

    CHECK_CUDA_ERRORS(cudaFreeAsync(d_receiver_pixel_pos, stream_));
    CHECK_CUDA_ERRORS(cudaFreeAsync(d_pixel_angles, stream_));
    if (!external_stream_)
    {
        CHECK_CUDA_ERRORS(cudaStreamSynchronize(stream_));
    }
}
